利用近期新冠疫情於台灣之確診者足跡製作熱點圖(核密度法,Kernel Density Estimation),並與台北市人口密度進行套疊,以觀察確診者足跡熱點空間分布與台北市人口密度的關係。
import pandas as pd
import geopandas as gpd
from pyproj import CRS
import requests
import geojson
import folium
# STEP1 加入統計資料
taipei_pop = pd.read_csv(r"E:\Py_NTU 341\Homework\Data Folder\110年4月台北市各里人口數戶數.csv",encoding='big5')
taipei_pop.head()
| 區域 | 村里 | 村里數_現有門牌 | 村里數_戶籍登記 | 鄰數_現有門牌 | 鄰數_戶籍登記 | 戶數 | 人口數-合計 | 人口數-男 | 人口數-女 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 松山區 | 莊敬里 | 1 | 1 | 27 | 27 | 2056 | 5279 | 2548 | 2731 |
| 1 | 松山區 | 東榮里 | 1 | 1 | 33 | 32 | 3071 | 8076 | 3773 | 4303 |
| 2 | 松山區 | 三民里 | 1 | 1 | 32 | 32 | 2875 | 6760 | 3086 | 3674 |
| 3 | 松山區 | 新益里 | 1 | 1 | 18 | 18 | 1843 | 4569 | 2174 | 2395 |
| 4 | 松山區 | 富錦里 | 1 | 1 | 23 | 23 | 2125 | 5181 | 2369 | 2812 |
# STEP2 新建一個TownVillage欄位,為區域和村里的文字合併(作為後續關聯使用)
taipei_pop["t.v"] = taipei_pop["區域"] + taipei_pop["村里"]
#taipei_pop.head()
pop = taipei_pop[["t.v","人口數-合計"] ]
pop = pop.rename(columns={'人口數-合計': 'pop202104'})
pop.head()
| t.v | pop202104 | |
|---|---|---|
| 0 | 松山區莊敬里 | 5279 |
| 1 | 松山區東榮里 | 8076 |
| 2 | 松山區三民里 | 6760 |
| 3 | 松山區新益里 | 4569 |
| 4 | 松山區富錦里 | 5181 |
STEP1 找到全台村里界圖 村里界圖(TWD97_121分帶)
STEP2 選出台北市之範圍
STEP3 計算村里面積(計算人口密度用)
STEP4 新建一個TownVillage欄位,為區域和村里的文字合併(作為後續關聯使用)
STEP5 計算人口密度
# STEP1 選出台北市之範圍
tw_village = r"E:\Py_NTU 341\Homework\Data Folder\VILLAGE_MOI_121_202104\VILLAGE_MOI_121_1100415.shp"
tw_village = gpd.read_file(tw_village,encoding = "UTF-8")
#tw_village.head()
# STEP2 選出台北市之範圍
tp_village = tw_village.loc[tw_village['COUNTYNAME']=='臺北市']
tp_village.head()
| VILLCODE | COUNTYNAME | TOWNNAME | VILLNAME | VILLENG | COUNTYID | COUNTYCODE | TOWNID | TOWNCODE | NOTE | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 6952 | 63000080031 | 臺北市 | 文山區 | 樟新里 | Zhangxin Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((306484.609 2763701.609, 306442.073 2... |
| 6953 | 63000080037 | 臺北市 | 文山區 | 老泉里 | Laoquan Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((307698.415 2763909.981, 307784.205 2... |
| 6954 | 63000080032 | 臺北市 | 文山區 | 樟腳里 | Zhangjiao Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((307451.461 2763912.590, 307424.695 2... |
| 6955 | 63000080041 | 臺北市 | 文山區 | 樟文里 | Zhangwen Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((306530.118 2763845.236, 306483.024 2... |
| 6956 | 63000080043 | 臺北市 | 文山區 | 樟樹里 | Zhangshu Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((306790.748 2763911.814, 306732.793 2... |
# STEP3 計算村里面積(計算人口密度用)
tp_village2 = tp_village.copy() #避免py:853: SettingWithCopyWarning
tp_village2['Area'] = tp_village.area/10000
tp_village2.head(2)
| VILLCODE | COUNTYNAME | TOWNNAME | VILLNAME | VILLENG | COUNTYID | COUNTYCODE | TOWNID | TOWNCODE | NOTE | geometry | Area | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 6952 | 63000080031 | 臺北市 | 文山區 | 樟新里 | Zhangxin Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((306484.609 2763701.609, 306442.073 2... | 19.511016 |
| 6953 | 63000080037 | 臺北市 | 文山區 | 老泉里 | Laoquan Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((307698.415 2763909.981, 307784.205 2... | 377.403620 |
#做 Town + village
tp_village3 = tp_village2.copy()
tp_village3["t.v"] = tp_village["TOWNNAME"] + tp_village["VILLNAME"]
tp_village3.head(2)
#village = tp_village[["COUNTYNAME","TOWNNAME","t.v","geometry","Area"]]
#village.head()
| VILLCODE | COUNTYNAME | TOWNNAME | VILLNAME | VILLENG | COUNTYID | COUNTYCODE | TOWNID | TOWNCODE | NOTE | geometry | Area | t.v | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 6952 | 63000080031 | 臺北市 | 文山區 | 樟新里 | Zhangxin Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((306484.609 2763701.609, 306442.073 2... | 19.511016 | 文山區樟新里 |
| 6953 | 63000080037 | 臺北市 | 文山區 | 老泉里 | Laoquan Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((307698.415 2763909.981, 307784.205 2... | 377.403620 | 文山區老泉里 |
# STEP4 新建一個TownVillage欄位,為區域和村里的文字合併(作為後續關聯使用)
tp_density = tp_village3.merge(pop,on="t.v") #.fillna(0) 補0
tp_density.head(2)
| VILLCODE | COUNTYNAME | TOWNNAME | VILLNAME | VILLENG | COUNTYID | COUNTYCODE | TOWNID | TOWNCODE | NOTE | geometry | Area | t.v | pop202104 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 63000080031 | 臺北市 | 文山區 | 樟新里 | Zhangxin Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((306484.609 2763701.609, 306442.073 2... | 19.511016 | 文山區樟新里 | 5662 |
| 1 | 63000080037 | 臺北市 | 文山區 | 老泉里 | Laoquan Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((307698.415 2763909.981, 307784.205 2... | 377.403620 | 文山區老泉里 | 912 |
# STEP5 計算人口密度
tp_density["density"] = tp_join["pop202104"] / tp_join["Area"]
tp_density.head(2)
| VILLCODE | COUNTYNAME | TOWNNAME | VILLNAME | VILLENG | COUNTYID | COUNTYCODE | TOWNID | TOWNCODE | NOTE | geometry | Area | t.v | pop202104 | density | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 63000080031 | 臺北市 | 文山區 | 樟新里 | Zhangxin Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((306484.609 2763701.609, 306442.073 2... | 19.511016 | 文山區樟新里 | 5662 | 290.195032 |
| 1 | 63000080037 | 臺北市 | 文山區 | 老泉里 | Laoquan Vil. | A | 63000 | A11 | 63000080 | None | POLYGON ((307698.415 2763909.981, 307784.205 2... | 377.403620 | 文山區老泉里 | 912 | 2.416511 |
#查看資料分布
plt = tp_density["density"].plot(figsize=(10,5))
plt.grid("on", axis="y")
plt.grid("on", axis="x")
plt.set_title("density")
plt.set_xlabel("VILLNAME")
Text(0.5, 0, 'VILLNAME')
#敘述性統計
tp_density["density"].describe()
count 455.000000 mean 284.789686 std 188.315154 min 0.545540 25% 141.283916 50% 280.230794 75% 403.193070 max 1112.840866 Name: density, dtype: float64
#展圖
tp_density.plot(figsize=(6,6),column='pop202104', scheme="natural_breaks",cmap='Greens', k=10) #legend=True
tp_density.plot(figsize=(6,6),column='density', scheme="natural_breaks",cmap='Greens', k=10) #legend=True
<AxesSubplot:>
covid19 = r"E:\Py_NTU 341\Homework\Data Folder\0516covid19_point.shp"
covid19 = gpd.read_file(covid19,encoding = "big-5")
#covid19.head(3)
covid19[["Name","geometry"]]
| Name | geometry | |
|---|---|---|
| 0 | 4/28 桃園市某國民小學 | POINT Z (121.28966 25.04654 0.00000) |
| 1 | 5/1Pilot in Cafe | POINT Z (121.20416 25.00696 0.00000) |
| 2 | 5/1胡同彭家老舖 新疆拉麵 | POINT Z (121.29158 25.04643 0.00000) |
| 3 | 5/1華泰名品城 | POINT Z (121.21327 25.01412 0.00000) |
| 4 | 5/3 台北福華大飯店 | POINT Z (121.54297 25.03747 0.00000) |
| ... | ... | ... |
| 272 | 教育部永續校園_新北市秀山國民小學 | POINT Z (121.52114 24.99516 0.00000) |
| 273 | 5/7 患者自行前往 診所就醫 | POINT Z (121.52060 25.00017 0.00000) |
| 274 | 三水街 | POINT Z (121.50136 25.03608 0.00000) |
| 275 | 西園路一段 | POINT Z (121.49945 25.03709 0.00000) |
| 276 | 雙和醫院-患者自行前往 | POINT Z (121.49381 24.99130 0.00000) |
277 rows × 2 columns
# Get x and y coordinates for each point
covid19["x"] = covid19["geometry"].apply(lambda geom: geom.x)
covid19["y"] = covid19["geometry"].apply(lambda geom: geom.y)
# Create a list of coordinate pairs
locations = list(zip(covid19["y"], covid19["x"]))
# 創建Choropleth map階級區分圖所需的Geo-id(每行必須具有唯一的標識符)
tp_density['geoid'] = tp_density.index.astype(str)
from folium.plugins import HeatMap
m = folium.Map(location=[25.04, 121.53], tiles = 'cartodbpositron', zoom_start=12, control_scale=True)
folium.Choropleth(geo_data = tp_density,
data = tp_density,
name='Density',
columns=['geoid','density'], #資料欄位
key_on='feature.id',
fill_color='Greens', #Choropleth map色塊
line_color='white', #Choropleth map線色
line_weight=0, #Choropleth map線粗
legend_name= 'Density').add_to(m) #圖例名稱
#工具提示->顯示每個功能的值
folium.features.GeoJson(tp_density,
name='Density label',
style_function=lambda x: {'color':'transparent','fillColor':'transparent','weight':0},
tooltip=folium.features.GeoJsonTooltip(
fields=['t.v','density'],
aliases = ['區域村里','人口密度(人/公頃)'],
labels=True,
sticky=True #顯示資訊方式不依樣(黏著度)
)
).add_to(m)
HeatMap(locations,name='HeatMap Covid-19 Case',radius=30).add_to(m) #假設熱點半徑為30km
# 加入圖層開關功能 layer control
folium.LayerControl().add_to(m)
m
outfp = "E:\Py_NTU 341\Homework\covid19heatmap&pop_density.html"
m.save(outfp)